## Analysis of Welfare Bills in subsection "Evidence of a Constituency Connection"
## Final version JBS 3 August 2017

rm(list=ls())

library(mvtnorm)
library(arm)

path <- "~/Dropbox (Personal)/Rebel Summaries/APSR_SKLLO_Repfiles/"

load("RawData/welfarerebels.Rda")

mod1 <- glm(rebel ~ bens00 , family=binomial(link='logit') , data=mps97)
summary(mod1)


sims<-1000

# Take random draws
dist<-rmvnorm(sims, coef(mod1), as.matrix(vcov(mod1)))

pp1<-numeric(sims)
pp2<-numeric(sims)

for (j in 1:sims){
  pp1[j]<-dist[j,1]+dist[j,2]*quantile(mps97$bens00 , 0.1)
  pp2[j]<-dist[j,1]+dist[j,2]*quantile(mps97$bens00 , 0.9)
 } 

pe1 <-invlogit(mean(pp1))
lo1 <-invlogit(quantile(pp1, 0.025))
hi1 <-invlogit(quantile(pp1, 0.975))

pe1

pe2 <-invlogit(mean(pp2))
lo2 <-invlogit(quantile(pp2, 0.025))
hi2 <-invlogit(quantile(pp2, 0.975))

pe2

fdiff <- mean(invlogit(pp2)-invlogit(pp1))
fdifflo <- quantile(invlogit(pp2)-invlogit(pp1), 0.025)
fdiffhi <- quantile(invlogit(pp2)-invlogit(pp1), 0.975)

fdiff